Exactly solving the Kitaev chain and generating Majorana-zero-modes out of noisy qubits

Majorana-zero-modes (MZMs) were predicted to exist as edge states of a physical system called the Kitaev chain. MZMs should host particles that are their own antiparticles and could be used as a basis for a qubit which is robust-to-noise. However, all attempts to prove their existence gave inconclusive results. Here, the Kitaev chain is exactly solved with a quantum computing methodology and properties of MZMs are probed by generating eigenstates of the Kitev Hamiltonian on 3 noisy qubits of a publicly available quantum computer. After an ontological elaboration I show that two eigenstates of the Kitaev Hamiltonian exhibit eight signatures attributed to MZMs. The results presented here are a most comprehensive set of validations of MZMs ever conducted in an actual physical system. Furthermore, the findings of this manuscript are easily reproducible for any user of publicly available quantum computers, solving another important problem of research with MZMs—the result reproducibility crisis.

In this Section I will briefly outline the procedure for obtaining the ground state of quadratic Hamiltonians relying heavily on steps described in Appendix A and Section IV of Ref. [1]. I will try to use the same notation and follow exactly the same steps. The reader is refereed to the original study for a more detailed discussion. This procedure is implemented in Google's Quantum AI Cirq in an exact numerical fashion [2]. Afterwards, The circuit obtained in Google's Quantum AI Cirq is converted to a Qiskit form and the code for that is given in 10.5281/zenodo.6323467.
The procedure is as follows: 1. Find A the Majorana representation of a general quadratic Hamiltonian H = i 2 f T Af . 2. Bring the matrix A into the Schur form and permute indices such that the upper block ε is diagonal with eigenvalues in increasing order perm(Schur(A)) = RAR T = 0 ε −ε 0 . (S1) 3. Obtain W = Ω † RΩ, where The form of W is in general 4. Find such unitary transformation U for which V W L U = [0 1], where V is a general unitary matrix, . 5. U can be decomposed into products of X gates on the last qubit and a combination of RY XXY and single qubit Z gates -and it is a representation of a unitary operator which diagonalizes a general quadratic Hamiltonian.

S2 -THE WINDING NUMBER -A TOPOLOGICAL INVARIANT FOR A FINITE KITAEV CHAIN
In this Section I will give an answer to the question of the boundary of the topological phase of finite Kiteav chains. To paraphrase it more specifically and in a more simple manner: when t = −∆ and n is finite, for which value of µ does a topological-non-topological transition occur in the Kitaev Hamiltonian?
The starting point is the derivation is n-site Kitaev Hamiltonian Here the capital K denotes the wavenumber while the lowercase k denotes the site in the Kitaev chain. Setting µ k = µ, t jk = t and ∆ jk = ∆ the following Hamiltonian is obtained One can draw a conclusion that the behavior of the infinite Kitaev chains in K-space is the same as those of finite chains with a modified tunneling and superconducting pairing t →t and ∆ →∆. When ∆ = −t the topological condition in infinite Kitaev chains is µ < 2t [4,5]. Given that topological properties are fully determined from properties of the K-space Hamiltonians [6] the topological condition for a finite Kitaev chain is The potentially topological state of matter in a finite Kitaev chain can be further described by a winding number ν [6]. Similarly with a Chern number the winding number is a topological invariant of a given system [5,6]. In the case of the Kitaev chain, ν = ±1 would indicate a topological phase and ν = 0 a topologically trivial phase. The winding number is defined as [5,6] where ∂ K w(K) is the winding number density and w(K) = Arg[2∆ sin (K) + i(µ + 2t cos (K))]. Here Arg stands for a non-principal value of the argument. Note that the winding number is a Z topological invariant and it was long thought that it cannot be used to characterize Kitaev chains. Work of Ref. [7] changed this and proved the equivalence of Eq. (S9) and a more adapted expression for Kitaev chains. The integral in Eq. (S9) is difficult to compute due to the fact that the non-principal value of the argument cannot be found in an analytical fashion. Still, a numerical routine can be developed, by unwinding an array of different values of w(K) in the region [−π, π]. In the results given at Fig. S1 we compare results of Eq. (S8) with a numerical solution of the integral in Eq. (S9) obtained by numerically unwinding w(K) in the region [−π, π] by Python's numpy unwind function.
In Fig. S1 we see that even the shortest possible Kitaev chain of only 2 sites, remains topological for µ < t = 1 and that as the length of the chain is increased the topological-non-topological boundary approaches that of infinite chains µ < t = 2 [4]. In the case of the experiment in the main body of the paper n = 3, µ cr = 4/3 and so the system remains topological in 75% of values of µ as compared to an infinite Kitaev chain. This means that the 9 points corresponding to low value of µ are in the topologically non-trivial regime in the main body of the paper.
Finally, two-qubit gates error rates of contemporary quantum computers (< 0.75%) and are the most dominant source of errors on current-day quantum devices. We can try to estimate are such constrains on error rates of two qubit gates preventing one to achieve an MZM state on a quantum computer.
To measure that we define a parameter which represents a difference between the two-qubit gate angles at µ = 0 + and µ ≈ 1.33 normalized to π for every two qubit gate in the circuit preparing the ground state of the MZMs (Fig. 1 -main body of the paper) (S10) In Tab I we see that circuit is remaining in the topological regime even if the uncertainty in generating the two qubit gate is on the order of ∼ 0.75% outlining the possibility of creating MZMs on modern gate-based quantum computers. The exception to this is the first two qubit gate between qubits q 1 and q 2 which has to have a value of ±π/2 for a broad range of µ. However, the two-qubit RY XXY gate has a special symmetry point at θ ij = π/2 in which it is equivalent to two single qubit gates Z-gates which have much lower error rates (on the order or < 0.1%)

S3 -A NOISE MODEL
Noise in a quantum system governed by the Hamiltonian H and described by a density matrix ρ is given by  δ absolute difference between the two-qubit gates angle of the initial topological and trivial regions. The circuit has the same form like in the main body of the paper in Fig. 1.
(S12) Here, Γ ij is the dissipation rate between the channel ij, and L † ij = |i j| and L ij = |j i| are Lindblad dissipation operators. In general if i = j the terms in the Lindblad equation describe relaxation processes happening on a characteristic relaxation timescale T 1 and if i = j these terms correspond to pure dephasing processes described by a characteristic time T * 2 . In general T * 2 T 1 -pure dephasing dominates over relaxation -and this is the assumption I adopt in the remainder of the Section.
The solution of a Lindblad equation in the presence of only pure dephasing is the diagonal part off the density matrix remains unchanged while the off-diagonal terms are decaying with a factor of exp (−τ /T * 2 ) (see Section "Phase damping" in Ref. [9]).
By repeating the procedure from the main body of the paper -obtaining the quantum circuit with OpenFermion and then getting the ground state wavefunction |ψ in a form of a numerical vector we can obtain the density matrix ρ(0) = |ψ ψ|. After obtaining the density matrix and multiplying its off-diagonal terms by exp (−τ /T * 2 ) we obtain the noisy density matrix of the systemρ (this simulates idle noise processes). The noisy expectation value of the Hamiltonian can be obtained as E = Tr(ρ H). (S14) Recent experimental and theoretical studies [10][11][12][13] indicate that the dephasing time of superconducting qubits varies in time. A minimal model in describing that is assigning random fluctuations of the dephasing time between measurements of different eigenstates of energy at different values of µ, T * 2 → T * 2 + ab where a = [−1, 1] is a random number and b T * 2 . In Fig. S2 I show what is the effect of pure dephasing on the spectrum of MZMs. As the duration of the quantum circuit is increased, the circuit dephases more (first to third subfigure) the topologically trivial states (E = ±2 at µ = 0 + ) move more towards zero in the BdG spectrum (black lines). Furthermore, variations in the dephasing time between the realization of the experiment yield a qualitatively similar spectrum to the one measured (green symbols in the 3rd subfigure vs. measured data -4th subfigure). This effect is visually similar to a slight violation of particle-hole symmetry. Noise in the quantum computing device could be understand as performing an MZM experiment with uncertainties in tunneling, chemical potential and superconducting pairing like that one would encounter in a realistic system.

S4 -THE MAJORANA EDGE CORRELATION FUNCTION
In the thermodynamic limit, the Majorana edge correlation function is expected to decay with a quadratic dependence in µ [14]. As comparing a n = 3 site Kitaev chain with a theory developed for an infinitely long chain is infeasible, I will just give mean square fits of the decay of the absolute value of the Majorana edge correlation function | iγ 1 γ k |, with a goal of outlining the fact that Majorana edge correlations do indeed decay as a function of µ. The polynomial fit to data for the first run (left) for state [0] is 0.57 − 0.18µ + 0.03µ 2 and for the state [1,2] is 0.56 − 0.18µ + 0.03µ 2 with respective residuals 0.03 and 0.02.
The polynomial fit to data for the second run (right, main body of the paper) for state [0] is 0.63 − 0.14µ + 0.03µ 2 and for state [1,2] 0.6 − 0.09µ + 0.01µ 2 with residuals of 0.01 and 0.03 respectively. The device was calibrated between runs.

S5 -COMPARISON BETWEEN RUNS OF HARDWARE
In this section I compare two different experiment runs on IBM Santiago in FIG. S3. Although the different executions vary slightly quantitatively depending on the calibration of the device all qualitative features remain present, especially the possible topological degeneracy. In Fig. S4 I display a comparison between a BdG spectrum and two different realizations of the experiment. Both experiments show excellent qualitative agreement with the BdG prediction. One can quantify this by defining a mean absolute error as ME = ( i |x i − y i |) /m, where m is the total number of measurements/predictions, x i is the ith measurement of energy from the quantum device and y i is the ith prediction of the BdG Hamiltonian. I find a ME = 0.127 for the left figure and ME = 0.129 for the right figure for −1.87 ≤ y i ≤ 1.87.

S6 -A 4-SITE KITAEV CHAIN
The formalism presented here allows the treatment of longer Kitaev chains with quantum simulators. Actual quantum computing preparation of Kitaev states should somewhat still be possible with currently available quantum computing hardware, however the author of this paper has no access to state of the art quantum computers and I will thus only focus on results from quantum simulators in the remainder of this Supplementary Material. In